# functions for difference in means (pre-test)

diff_mean <- function(data, treat, questions){
  YT_ls <- get_YT(data=data, treat=treat, questions=questions)
  if(!is.null(YT_ls$Y_intl)){
    intl <- diff_mean_CI_helper(Y=YT_ls$Y_intl, treat=YT_ls$T_intl)
  }else{
    intl <- NULL
  }

  if(!is.null(YT_ls$Y_const)){
    const <- diff_mean_CI_helper(YT_ls$Y_const, YT_ls$T_const)
  }else{
    const <- NULL
  }

  if(!is.null(YT_ls$Y_ngo)){
    ngo <- diff_mean_CI_helper(YT_ls$Y_ngo, YT_ls$T_ngo)
  }else{
    ngo <- NULL
  }
  out <-rbind(intl, const, ngo)
  colnames(out) <- c("lower", "point", "upper") 
  return(as.data.frame(out))
}



get_YT <- function(data, treat, questions){

  Y <- apply(as.matrix(data[,questions]), 2, as.numeric)
  Y <- apply(Y, 1, sum, na.rm =T)
  Y <- ifelse(Y == 0 | Y == 99, NA, Y)

  Y_intl <- Y[treat$control == 1 | treat$intl == 1]
  T_intl <- treat$intl[treat$control ==1 | treat$intl == 1]
  
  Y_const <- Y[treat$control == 1 | treat$const == 1]
  T_const <- treat$const[treat$control ==1 | treat$const == 1]

  Y_ngo <- Y[treat$control == 1 | treat$ngo == 1]
  T_ngo <- treat$ngo[treat$control ==1 | treat$ngo == 1]

  return(list(Y_intl=Y_intl, T_intl=T_intl,
              Y_const=Y_const, T_const=T_const,
              Y_ngo=Y_ngo, T_ngo=T_ngo))
}

diff_mean_helper <- function(Y, treat){
  # obtain different in means estimator and related statistics 

  Y_trt <- Y[treat == 1]
  Y_ctr <- Y[treat == 0]
  tau_hat <- mean(Y_trt, na.rm=T) - mean(Y_ctr, na.rm=T)
  n_trt <- sum(!is.na(Y_trt))
  n_ctr <- sum(!is.na(Y_ctr))
  tau_se <- sqrt(var(Y_trt, na.rm=T)/n_trt + var(Y_ctr, na.rm=T)/n_ctr)
  z_val <- tau_hat/tau_se
  p_val <- (1-pnorm(abs(z_val)))*2
  CI <- c(tau_hat + qnorm(0.025)*tau_se, tau_hat, tau_hat + qnorm(0.975)*tau_se)
  return(list(tau_hat=tau_hat, tau_se=tau_se, z_val=z_val, p_val=p_val, CI=CI))
}

diff_mean_CI_helper <- function(Y, treat){
  Y_trt <- Y[treat == 1]
  Y_ctr <- Y[treat == 0]
  tau_hat <- mean(Y_trt, na.rm=T) - mean(Y_ctr, na.rm=T)
  n_trt <- sum(!is.na(Y_trt))
  n_ctr <- sum(!is.na(Y_ctr))
  tau_se <- sqrt(var(Y_trt, na.rm=T)/n_trt + var(Y_ctr, na.rm=T)/n_ctr)
  CI <- c(tau_hat + qnorm(0.025)*tau_se, tau_hat, tau_hat + qnorm(0.975)*tau_se)
  return(CI)
}


